Exact dynamics of dissipative electronic systems and quantum transport: 
Hierarchical equations of motion approach 
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A quantum dissipation theory is formulated in terms of hierarchically coupled equations of mo- 
tion for an arbitrary electronic system coupled with grand canonical Fermion bath ensembles. The 
theoretical construction starts with the second-quantization influence functional in path integral 
formalism, in which the Fermion creation and annihilation operators are represented by Grassmann 
variables. Time-derivatives on influence functionals are then performed in a hierarchical manner, on 
the basis of calculus-on-path-integral algorithm. Both the multiple-frequency-dispersion and the 
non-Markovian reservoir parametrization schemes are considered for the desired hierarchy construc- 
tion. The resulting formalism is in principle exact, applicable to interacting systems, with arbitrary 
time-dependent external fields. It renders an exact tool to evaluate various transient and stationary 
quantum transport properties of many-electron systems. At the second-tier truncation level the 
present theory recovers the real-time diagrammatic formalism developed by Schon and coworkers. 
For a single-particle system, the hierarchical formalism terminates at the second tier exactly, and 
the Landuer-Biittiker's transport current expression is readily recovered. 

PACS numbers: 72.10.Bg, 05.30.-d 



I. INTRODUCTION 

The aim of this work is to establish an exact quan- 
tum dissipation theory (QDT), with which various quan- 
tum transport properties could in principle be evaluated 
without approximations. Quantum transport through 
nanosystems has conventionally been studied via the 
Landauer-Buttiker scattering theory^ and nonequilib- 
rium Green's function (NGF) formalism. 2 These ap- 
proaches are however basically single-particle theo- 
ries. The fundamental physics of quantum transport 
through an interacting system has been understood 
based on methods applicable to certain limits; for ex- 
ample, in cither the weak or strong Coulomb interac- 
tion regime Schon and coworkers have devel- 
oped a formulation based on real-time diagrammatic 
technique^ ' 10 ! 11 ' 12 This approach could in principle 
be exact; however, to identify all diagrams is itself 
formidable, and there are practically no feasible ways to 
overcome this difficulty for a general non-Markovian sys- 
tem in the presence of arbitrary Coulomb interaction. 
Other promising nonpcrturbative methods, such as the 
numerical renormalization group approach ) ' ' are yet 
to be extended to dynamical quantum systems. 

Quantum transport is studied with the situation, 
where a "device" , such as a semiconductor quantum dot 
or organic molecule under investigation, is connected to 
electrodes under applied bias voltages. This situation 
can be well described in the framework of QDT. The 
latter concerns the fundamental formulation that gov- 
erns the dynamics of a quantum open system. The pri- 
mary quantity in QDT is the reduced system density op- 
erator, p(t) = tr B pr(*)- Here, pr(t) denotes the total 
density operator of system-bath composite; tr B the par- 
tial trace over electrode bath degrees of freedom. Quan- 
tum transport based on QDT approach has been formu- 



lated extensivelyr^iiiii^ii^i^iSii 2 ^! 2 ^ including the afore- 
mentioned real-time diagrammatic formalismi&^ i 10 ' 11 ' 12 
This approach has the advantage of its generality, since 
different scattering processes can be handled in a unified 
manner, and transient dynamics can be studied readily. 
However, the QDTs used in quantum transport by far are 
all perturbative in nature; most of them are only of the 
second order in the system-bath coupling. 18 ' 19 ' 20 ' 21 ' 22 i 23 
A second-order theory is also related to the sequential 
tunneling regime. Other limitations include the moder- 
ately high temperature/ voltage and quasi- broadband (or 
quasi-Markovian) approximation. Quantum transport 
with cotunneling processe s 24 ' 25 has been studied with 
the fourth-order QDTs, constructed via the standard 
projection operator technique ! 26 ' 2 ' 28 A self-consistent 
Born approximation to nonperturbative QDT has also 
been proposed to recover such as the nonequilibrium 
Kondo effect in a model system* 2 ^, Approximations in- 
volved in the existing quantum transport theories, on 
the basis of either QDT or NGF formalism, are subject 
to ever increasing challenge, due to especially the emerg- 
ing fields of quantum measurement and quantum infor- 
mation processes . 29 ' 30 ! 31 ' 32 

This work continues our recent effort on the devel- 
opment of QDT formalism i 33 ' 34 ' 35 ' 36 The most relevant 
one is Ref. [36|, in which a nonperturbative theory was 
constructed for open quantum system, interacting with 
Boson-like grand canonical bath ensembles. The present 
paper exploits the second-quantization field theory that 
properly treats fermionic transfer coupling processes. 

The remainder of the paper is organized as follows. In 
Sec. UH we specify the electron transfer coupling Hamilto- 
nian and the fluctuation-dissipation theorem to be used 
in the later development. In Sec. lIII Al we revisit the 
influence functional path-integral formalism, with the 
Fermion field representation that involves Grassmann 
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variables. The derivation is detailed in AppendixfX] Pre- 
sented in Sec. lIIIBl is the QDT-based expression of tran- 
sient transport current. In Sec. HV| we first consider a 
differential form of QDT, by exploiting the calculus-on- 
path-integral (COPI) algorithm^^^^L^^ We then 
develop a multi-frequency-dispersed hierarchical equa- 
tions of motion (MFD-HEOM) formalism of QDT. The 
final MFD-HEOM results are summarized in SecfV] The 
present theory is in principle exact, applicable to arbi- 
trary quantum transport systems. It is shown in Sec. lVIl 
to recover the celebrated Landauer-Biittiker's transport 
current expression for single-particle systems. The in- 
volving reduced single-particle density matrix dynam- 
ics is derived in Appendix |B] Moreover, it is demon- 
strated in Sec. lVIII that at the second-tier truncation the 
present MFD-HEOM theory recovers the real-time dia- 
grammatic formalism developed by Schon and cowork- 
ers for Coulomb interaction systems. In Sec. I Villi we 
present an alternative but equivalent HEOM formalism 
via parametrization, with the derivation detailed in Ap- 
pendix [U] Finally, we conclude this work in Sec. lIXl 



II. STOCHASTIC COUPLING HAMILTONIAN 
AND QUANTUM STATISTICAL MECHANICS 

A. Stochastic transfer coupling Hamiltonian 

Consider an electron (or spin) transport setup, in 
which a multi-level system such as a semiconductor quan- 
tum dot or molecule is in contact with electrodes (labeled 
by index a). Each electrode serves as an electron reser- 
voir and is treated as a grand canonical Fermion bath 
ensemble. The total system-electrodes composite Hamil- 
tonian assumes 



(2.1a) 



with the noninteracting electrons in the a-electrode, 

h a = y^ y e ak dl ik d a k, (2-lb) 
fe 

and the system-electrode transfer coupling, 

H ' a = J2 W rf lfe«M + H.c. (2.1c) 

H.c. stands for the Hermitian conjugate; aj^ and d ak (a^ 
and dak) ar e the creation (annihilation) operators for an 
electron (or spin) in the specified spin-orbital of the sys- 
tem and the a-electrode, respectively. The system Hamil- 
tonian, H = H(t; {ajj, a^}) in the right-hand-side (rhs) of 
Eq. (|2. la|) . is rather arbitrary, including Coulomb inter- 
action and coupling with external time-dependent fields 
such as the electric field induced by a pulsed laser. 

Throughout this work, we set h = 1 and f3 a = 
l/(k B T a ), with k B the Boltzmann constant and T a the 



temperature of a-electrode. Electrons in bare electrode 
in steady state, either equilibrium or noncquilibrium, are 
described by the grand canonical density operator, 



Pi 



-0 a (h a -n a N a ) 



tl'r 



(2.2) 



Here, fj, a denotes the chemical potential of the a- 
electrode in steady state; N a = d ak d a k is the particle 
number operator of electrons in a-electrode. It satisfies 
[N a ,h B ] = 0, where h B = ^ Q h a . Denote also 



(0> B EEtr B (dp B ); 



ith 



IL 



(2.3) 



To describe the stochastic nature of the transfer cou- 
pling, consider Eq. (|2.1c| in the reservoir /i B -interaction 
picture: 



with 



fUt) 



k 



-ihnt 



(2.4) 



(2.5) 



being the stochastic interaction bath operators. They 
satisfy the Gaussian statistics with the Wick's theorem 
for thermodynamic average. Also note that (f afJ- (t)) B = 

(4(')tW>B = 0, and </t M (i)jU(r)} B = if a ± 
a' . As results, the effects of reservoirs on the reduced 
system can be completely determined by the two-time 
correlation functions, 

C+ fiV (t-r) = (fUt)faAr)) B , (2.6a) 
C-^(t-r) = (U(t)fl(r)) B . (2.6b) 

It follows immediately the time-reversal symmetry and 
the detailed-balance relations^ 

[Ct^(t)V = Ctpi-t) = e ±/3 "^C^(t - i(3a). (2.7) 

Physically, C+ M „ (t) describes the processes of electron 
tunneling from the reservoir a into the system, while 
C~^ v {t) describes the reverse events. Apparently, the 
reservoir correlation functions are diagonal with respect 
to spin indices; i.e., = 0, if [i and v belong to 

different spins. By far, the reservoir states are assumed 
to be time-independent. The resulting correlation func- 
tions defined in Eq. (|2.6p satisfy the stationary condition 
of C° v {t,r) = C° v (t — t). Nonstationary correlations 
for the case of time-dependent chemical potentials ap- 
plied on electrodes will be considered in Sec. Ill Cl see 



B. Fluctuation— dissipation theorem 

For sake of bookkeeping, we shall hereafter use also a 
to label the + or — , while a = —a the opposite sign of 
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a. Thus, Eq. (|2.T[) can be recast as 

[C^MY = Cau^-t) = e a ^C^ v (t - ip a ). (2.8) 

Denote also for either the creation (a+ = ajj or the 
annihilation (eC = a M ) operator for the specified spin- 
orbital state of system. 

Introduce now the spectrum functions r^ MJ/ (w) via 



C^ v (t)= I dwe^'r^H. (2.9) 



The above definition is consistent with the fact that 
the correlation functions defined in Eq. (|2.6[) are of 
= and C-^it) = i£>„,(i), in rela- 

tion to the self-energy functions in the Green's function 
technique;- The present index scheme is however more 
convenient in the construction of QDT formalism. The 
frequency-domain counterparts of Eq. (|2.8[) are 



rVM = [rs,»]* = e-^-^r*^). (2.10) 

The spectrum functions are of positivit y 33 ' 36 satisfying 
> and Tl^{w)T° avv {u) > |r^»| 2 . 
To express the detailed-balance relation in terms of 
fluctuation-dissipation theorem (FDT), let us consider 
the interaction spectral density functions, defined as 

1 C°° 

J apv (u)) = - / dte^-^{{f a ,{t)JUr)}) u . 

(2.11) 

For the present model of linear coupling with noninter- 
acting bath reservoir, it can be evaluated as J a ^v{^) — 
Ektlk^akJi" ~ fat)- Together with Eqs. jOJ, (H}, 
and the first identity in Eq. (|2.10p , we have 

= + r+ v/t M = -C M (w). (2.12) 

The detailed-balance relation [i.e. the second identity in 
Eq. (f2~T0]) ] leads then to 



r i^M = /a( w )-W(w), 



(2.13a) 

= /a M J <wM, (2.13b) 



with /+(w) = 1 — f a (ui) = f a (u)) being the Fermi distri- 
bution function; i.e., 



1 



1 + e <r ft*( a ' _ ' 1 <*) ' 



(2.14) 



By settin g J+^ u = J aUfJ , and J QAI „ = J a „ M , one can write 
Eq. (l2~T3l) as 



dco 



'JIM 



1 + e o'/3 Q (w-^ ) ' 



(2.15) 



This is the FDT in the fermionic grand canonical ensem- 
bles. It relates the correlation functions to the spectral 
densities. 



C. Correlation functions in the presence of 
time-dependent chemical potentials 

We shall also be interested in the transient dynamics 
of charge transport under time-dependent bias voltage. 
Its effect can be described by rigid homogeneous time- 
dependent shifts of the conduction bands of electrodes; 
i-e., e akfl (t) = e akfl + A a (t) and fi a (t) = p a + A a (t), 
so that the occupation on each state is unchanged^ As 
results, the nonstationary correlation functions are 



exp 



dt'A a (t') 



(2.16) 



for t > t. Note that A Q (i) in this work represents time- 
dependent chemical potential, on top of the constant part 
(fi a — ^° q ) applied on a-electrode. 



III. QUANTUM TRANSPORT VERSUS 
DISSIPATION 



A. Influence functional in path integral formalism 

We shall be interested in a QDT-based transport for- 
mulation. Let us first revisit the path-integral (PI) in- 
fluence functional expression. It serves as the start- 
ing point for the development of HEOM formalism of 
The quantity of interest here is the 
reduced density operator, defined as the trace of the to- 
tal density operator over the bath subspace; i.e., p(t) = 
tT B pr(i). The quantum dissipation starts with the initial 
condition that the system and bath were initially uncor- 
rected: 



pt(*o) = p(to)p° B ; 



tc, 



(3.1) 



This initial factorization ansatz does not cause approx- 
imation, as long as the initial time is set to infinite 
past i 33 ' 40 ' 41 We will come back to this issue in Sec. lIXI 

The key variation from our previous work with Boson- 
like reservoir— is the second quantization that leads to 
the Grassmann variables for Fermion fields rather than 
the ordinary c-numbers in the PI representation. 42 ! 43 To 
proceed, let {IV 1 )} be a second-quantization basis set in 
the system subspace, and ip = (ip, if)') for short, so that 
p(ip,t) = p(ip,ip',t). Denote U(t,to) as the reduced 
Liouville-space propagator, by which 

p{t)=U(t,t )p(t ). (3.2a) 

Its PI expression in the ^-representation reads^ 4 - 

rV»[t] 

'•<Mtoj 



/•V>[t] 

U{ip,t;ipo,t ) = Vipe iS[ ^T^}e- lS ^ ] . (3.2b) 



S[ip] is the classical action functional of the reduced sys- 
tem, evaluated along a path ip(r), with the constraints 
that two ending points ip(to) — ipo and if)(t) = if) are 
fixed. 
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The key quantity in PI expression is the influence func- 
tional, T in Eq. (|3.2b|) . For the present model of lin- 
ear coupling with noninteracting electron reservoir, it 
can be formally evaluated by using the Wick's theorem 
for the thermodynamic Gaussian average, implemented 
with Grassmann algebra i 42 ^ 43 The final results read [cf. 
Eq. lf£6]) ] 



^]=exp|-^ dr^[r;{^}]| 



(3.3a) 



with (denoting et+ = and 



fJ' fj, "-,u — VflJ 



K[t;{i>}} ee* ]T Al[i>(t)]B^(t; {!/>}). (3.3b) 



a,/j.,cr 



Here, A° and SJ„ are Grassmann variables, denned as 



>«(t)] = a>(t)]+a^(t)], 



(3.4) 



and 



B^Ct; {</>}) = -i[BS„(*5 {V}) ~ B^(t; {V/})], (3.5) 
with 

^(t; W) = E / drC^TKMr)}, (3.6a) 

„ "'to 

B^(t;{V'}) = E fdrCl; v {t,r)alW{r)]. (3.6b) 

„ "'to 

Note that B'^ = (B*Jt leads to B^ = (B^)t. This 
property corresponds to the Hermitian conjugate relation 
of auxiliary density operators that will be discussed later; 
see Eqs. (|3.12[) . (|5.5p or (|8.8p . The interaction reservoir 
correlation functions in Eq. (|3.6p can be nonstationary 
to support the study of transient current under time- 
dependent bias voltage applied to electrodes. 

All variables involved in Eqs. (|3.4[) and (|3.5p stem 
from second-quantization of fermion operators. They are 
Grassmann variables in PI formalism, 42 following the an- 
ticommutation relation, such as B^a^^ = — a a u B a a ll . Ap- 
parently, the dissipation functional 1Z [Eq. (|3.3b|) ]. as it 
consists of bi-fermionic variables, and the influence func- 
tional T [Eq. (|3.3ap ] remain as ordinary c-numbers. 

The time derivative on T reads 

d t T = -TIT = -iJ2 -KK^ = ~i E-^X- ( 3 - 7 ) 



Introduce here are a set of first-tier auxiliary influence 
junctionals (AIFs), 



■pa _ ger j: 



(3.8) 



They are Grassmann variables, for which [cf. Eq. (|3.4|) ] 

aIt° = <m)\Ku - ( 3 - 9 ) 



The identity of a*[i/>' (t)]^ = -^o*[^'(t)] for Grass- 
mann variables is used here. 

The first-tier auxiliary density operators (ADOs) can 
now be defined via [cf. Eq. (|3.2p ] 



(3.10a) 



with 



i>o[to] 



VtP e iS Mf°[tp]e- iS W' 



(3.10b) 

We can then recast Eq. (|3.7[) in terms of the reduced den- 
sity operator and its auxiliary ones as 



m = -i£p®-i^K,Pi ll (t)] 



(3.11) 



The first term in which Cp = [H, p] arises from the 
time derivative of classical action exponential term of 
Eq. (|3.2[) . The commutator in the second term can be 
recast as [a a ,Paa] = •AZPau tnat i s the operator-level 
counterpart of Eq. (|3.9[) . The Hermitian conjugate rela- 
tion implied in Eq. p.5p reads 



Po 



(3.12) 



This identity can also be verified via the equivalent 
operator- level definition of the first-tier ADO p^; see 
Ea. flgTB . 



B. Quantum transport current via the auxiliary 
reduced density operator dynamics 

In this subsection, we shall show that the transport 
current is directly related to the first-tier ADOs. The 
final expression for the transient current from a-lead to 
system reads (setting e = 1) 

/aW^^tr^+^K-a^W}. (3.13) 



Here, p^^ — yPau) are the aforementioned first-tier 
ADOs. The total current passing through the system 
from L to R electrode is then I(t) = Ih(t) — ^r(*)- 

We are now in the position to prove the current ex- 
pression, Eq. (|3.13[) . Let us start with the definition, 

I a (t) = -j t (N a ) = i([N a ,H T {t)]) T . (3.14) 

Here, ( • ) T = Tr T [ • p T (t)], with Tr T ee tr s tr B ; N a is the 
electron number operator in the electrode a, satisfying 
[N a ,H(t)} = [N a ,h B ] = 0. In Eq. (EUD, H T = H + 
H'(t) is the total composite Hamiltonian given in the 
/i B -interaction picture. One can readily obtain 

<[^(*K-°^W])t- ( 3 - 15 ) 



I a (t) 
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The equivalence between Eq. (|3 . 14[) and Eq. (|3.15p can 
be established by recognizing the fact that the first-tier 
ADOs defined in the previous subsection can be recast 
as 

^W=tr B [/^(t)p T W]. (3.16) 

In fact pr(t) = — i[ifr(*), Pr(t)] m this work reads [cf. 
Eq. (EH)] 

Applying now the definition of p(t) = tr B /3x(i), the trace 
cyclic invariance, and Eq. (|3.11[) . it concludes immedi- 
ately that Eq. (|3.16[) does amount to the first-tier ADO 
defined earlier. This equivalence can also be proved by 
the explicit PI expression of Eq. (|3.16p . 

We have therefore established the transport current 
expression [Eq. (|3. 13(1 ] . in terms of the reduced density 
operator, more precisely the first-tier ADO dynamics. 
Note that the charge current in Eq. (|3 . 1 3[) is of the form 
of I a (t) = Ylu^cttJ,(t)- It consists of the contributions 
from individual spin-orbitals that couple directly with 
the lead. Physically, it is useful to study individual con- 
tributions, for example spin-dependent current. 



IV. HIERARCHICAL EQUATIONS OF MOTION 
APPROACH TO QUANTUM DISSIPATION 

To complete the formulation, we shall develop hier- 
archically coupled EOM for evaluating both p and its 
ADOs, in the presence of arbitrary non-Markovian dis- 
sipation. Both MFD and parametrization schemes of 
hierarchy will be considered. The desired theory will 
be constructed via the calculus-on-path-integral (COPI) 
algorithmj 34 i 35 i 36 i 37 i 38 i 39 For clarity, we present in this 
section the derivation of MFD-HEOM formalism, while 
summarize its final results in the next section. The 
parametrization based HEOM formalism will be post- 
poned to Sec. lVIIll after its MFD equivalence is scruti- 
nized in relation with the Landauer-Buttiker's and the 
real-time diagrammatic formulations of quantum trans- 
port. 



A. Quantum dissipation theory via calculus on 
path integrals 

In the following development we omit the explicit PI 
variables dependence whenever it does not cause confu- 
sion. Adopt further the abbreviation j — {/ic} or {apcr} 
whenever applicable. So that Aj = AZ, Bj = B a ail , and 
Eq. (|3.7p is written as 

d t T= -i^A-jBjT. (4.1) 

3 



Note that Aj [Eq. (|3.4[) ] depends on the field at the fixed 
ending point ip(t) = ip of the path. The difficulty of 
the PI formalism is the evaluation of Bj [Eq. (|3.5p ] that 
involves at all to < r < t. To resolve the memory 
contained in Bj, consider its time derivative, 

d t Bj=Bj-iCj. (4.2) 

Here. 

Bj=B^ = -i[B^-B' a l], (4.3a) 
that is similar to Bj, but with [cf. Eqs. (|3.5[) and (|3.6[) ] 

££ M (*;M) = £ I drC^ v {t,T)al[i,{r)], (4.3b) 

and [denoting C»£ v = C^(0) = 0^(0); cf. Eq. (22)] 

Cj = ]T {cXMit)] cX<W(m ■ (4-4) 

Note that Cg^t, t) = C£ M „(0) used in Eq. g3]) is valid 
for either stationary or nonstationary reservoir correla- 
tion functions. Apparently, both Bj and Cj are Grass- 
mann variables. The former contains memory, while the 
latter depends only on the fixed ending point of the path. 

To continue the hierarchy construction, consider now 
the time-derivative on the first-tier AIF of Eq. (|3.8[) . Us- 
ing Eqs. (HH) and (@T2|), we have (noting that Tj = T a ap ) 

d t Tj = d t {BjT) 

= {Bj-%Cj)T-iY J KjAf®i>? 
j' 

= BjT - iCjT - i AyByBjT 

f 

J\ iC,:r ^.l,,/;,. (4.5) 

i' 

Note that Tj = BjT , which is a Grassmann variable, be- 
longs to a different class of first-tier AIF. We shall return 
to this issue in the next subsection. 

The last term in Eq. (|4.5p defines also the second-tier 
AIFs; i.e., 

r u b,b,t r !r (4.6) 

The second identity that implies Tjj = arises from 
the Grassmann anticommutation relation of Bj'Bj — 
-BjBj,. 

The second-tier AIFs are ordinary c-numbers. Thus, 
AyTjji in the last term of Eq. (|4. 5|) amounts to [noting 
that Aj, = A^, of Eq. (j3~4)) ] 

AyTjj, = a^im^jf + Fir°%>W(M- (4-7) 

It is in contrast to Eq. (|3.9p . which involves a first-tier 
AIF that is a Grassmann variable. 
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We are now in the position to summarize the above 
hierarchy construction. An n th -tier AIFs can be defined 
in general as 

-']" , (4-8) 

Note that exchanging any two indexes in a given AIF 
companies with a "— " sign. This leads to the property 
that the n indexes in J~j 1 j 2 ---j n should all be different; 
otherwise, the AIF would be zero. Denote also 

= (B jn B jn _, ■ -B n + ■■■ + B jn - ■ -B h % ) T. (4.9) 

Equations (|4.1|) and (|4.2|1 lead then to 

dtjj n) =^ n) ~i{C jn B jn _^ ••%+•■ ■ 

+ B, ■ ■ -B n C n )T - iJ^B jn - ■ ■B n AjBjT 

3 

n 

= ^)_ i J2 H n-kc jkB . n ,.. B . h+i 
k=l 

(4.10) 

The sum ^ runs over all j ^ jk', k = 1, • • • , n. Denote 
= B jn - ■ -B, . B,. ■ ■ -B h T, (4.11a) 

and 

•'•j.;' 1 B,B S . ■■■B s: :F. (4.11b) 
Equation (|4.1Q[) can now be recast as 

n 
fe=l 

-EVfi^ (4^) 

3 

Apparently, an odd-tier AIF is a Grassmann variable, 
while an even-tier one an ordinary c-number. This is the 
Grassmann parity that will affect the explicit expressions 
of the A- and C-terms in Eq. (|4.12[) in PI repression, and 
consequently their ADO counterparts in operator level; 
see Eqs. (l4~T5|) . 

B. Auxiliary density operators dynamics 

The n th -order ADO p ™ can now be defined via the 
AIF T\ n) as [cf. Eq. rfXTUD ] 

p^(t)=uj n \t,t )p(t ), (4.13a) 



with 

Jil>o[to] 

(4.13b) 

Similarly, the auxiliary operators p\ n \ p\ n , and 
are defined via , •^j'™ i an d respectively. 
The leading term in the n th -tier AIF such as or 
is of the (2n) th -order in the system-bath coupling; 
so is that in each individual p\ n ^ or pj . 

The EOM for p( n > can then be readily obtained via 
Eq. (|4~T2| and (|4~T3|) . Together with Eq. (|3~TTj) and set- 
ting p' -* = p, we have 

^-iCp-i^.AjpW, (4.14a) 

3 

n 

P f = -up™ +p\ n) -iEH^^pSr 1 ' 

/c=i 

-*EVl n+1) ; n>0 - ( 4 - 14b ) 

i 

Here, .Aj = AZ and C 3 = are superoperators. Their 
PI expressions that follow the Grassmann algebra are 
given in Eqs. (|3.4p and (|4.4p . respectively. As the Grass- 
mann parity associated with AIFs, the actions of these 
superoperators in Eq. (|4. 14[) are given as 

A- jP ^ = alp^ + {-) m p^a% (4.15a) 

and 

c 3 p (m) = E K» (m) - (-rcXp (m) <)- 

(4.15b) 

The n indexes in p\ n ' = py\ ■ should all be dis- 
tinct, due to the Grassmann or Fcrmion anticommuta- 
tion relation associated with the index swap. Recall that 
j = {a/icr}. As results, p {n) = 0, if n > 2N a N c . Here, 
N a denotes the number of electrodes; N c is the number 
of spin-orbitals that couple directly to the electrodes, [i.e, 
those of C a ^^(t, t) 7^ 0]. The factor 2 arises from the two 
possible signs of a = + and — . 

The EOM formalism presented in Eq. (|4.14p is finite 

but yet to be closed due to the fact that {/5j }, as defined 

via jf^ [Eq. IH35J)], involves both the B- and B-types 
functionals. The desired hierarchy should be constructed 
only via the £>-type functionals. To obtain a closed 
HEOM formalism, it requires that £>-type functionals 
[Eq. (O] be expressed in terms of the B's [Eq. (|33)) ]. 
This can be achieved via the extended Meier- Tannor pa- 
rameterization method ) 33 i 35 ' 36 i 41 ' 45 in which correlation 
functions C°v(t) are expanded in an exponential se- 
ries. This method was adopted in our previous develop- 
ment of QDT formalism for bosonic bath cases, in which 
the path integral variables are all c-numbersi 35 ' 36 We 
shall treat the parametrization approach to HEOM for 
Fermion reservoirs later; see Sec. I Villi and Appendix [Cl 
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C. The multiple-frequency dispersed hierarchy 
construction 

Proposed here is the MFD hierarchy scheme, in which 



where [noting that Bj k = Bj h {iOk,t) {ip})] 
^ (n) (a;;t;W) = 4v44^. 



(4.25) 



The proper MFD- ADO, <j>[ n) in Eq. (|4~T7)l . for closing the 



duj r . 



fijj'-jn ' ') ^n', t) ■ (4.16) hierarchy is now determined as 



The desired MFD-HEOM [cf. Eq. ([57EJ)] will be con- 
structed for a proper 



(«), 



(4.26a) 



with 



(4.17) 



Let us start with the time-independent chemical po- 
tential case, in which the reservoir correlation functions 
is stationary, C^^i, r) = C^^it — r), and can be ex- 
pressed in terms of spectrum functions as Eq. (|2.9p . In 
this case, Eq. (|3.6p can be recast as 



(4.18) 



(n) 



V>[*] 



, ! V^e iS ^Pi n) e- iS ^'\ (4.26b) 



with 



tf>(r)]. 



Note that the Grassmann anticommutation relation is 
now manifested not only in the j-indexes, but also in the 
associated frequencies; cf. Eq. (|4.25|) . 

The time-derivative on Eq. (|4.25|) can be obtained as 

n 
k=l 

-% J dwj2^" +1) (^^w)> ( 4 - 27 ) 



(4.19) 



The frequency-dispersed B'- is si milar to Bj, except w ith = J2k a k[^k + A OJk (t)]. The MFD-HEOM for 
that the T° oh the rhs of Eq. gH ijL£^««d by An) then b dn obtained- 
r a^ = j cf - the first identity of Eq. (123U1) ]. As re- ^ J J 
suits, Eq. (|3 !5[) can be expressed as 

V. HIERARCHICAL EQUATIONS OF MOTION 
WITH MULTIPLE FREQUENCY DISPERSION 

The final MFD-HEOM formalism at operator level 
reads as (setting <^ _1 ) = O (0) = and </> (0) = p) 



Bj = B^(t; {i>}) = J dioBj (u>, t; {</>}), (4.20) 

with Bj = -i{B 3 - B'j). 

In the presence of time-dependent chemical potentials 
applied on electrodes, the correlation functions include 
also the nonstationary phase factors [cf. Eq. (|2.16p ]. The 
frequency dispersed Bj -functional reads now 



. e tow(t-r) exp 



ai \ dt'A a (t') 



Bj = -i V f dr> 

x {r^>K[V>(r)] - K v ^)alW( T )]} . (4.21) 
It satisfies 

d t Bj = ia[uj + A a (t)]Bj + Cj(w), (4.22) 
with Cj(uj) = Cj(uj; {ip[i\}) being 

c i = E {K^Kim} - r*„>Khf'(t)]} • (4.23) 

V 

The n th -tier AIFs, T\ n) of Eq. gH), can now be dis- 
persed as 



n 

-^(-)»-% fc K)^- 1) K,i) 
fe=i 

-i^dw ^ ^^(w.o;,*). (5.1) 



Here, = {u>i,- ■ u k -i, u>k+ii ■ ;U n }, 



fc=l 



(4.24) 



and [cf. Eqs. and (1051) ] 

1/ 

- (-) m rt,>)* (m) <] 



(5.2) 



(5.3a) 



(5.3b) 
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The initial conditions to Eq. (|5.1[) are in fact those of 
steady-state solutions; see Sec. lIXI for details. The tran- 
sient current [Eq. (|3. 13|l ] is evaluated via the first-tier 
frequency-dispersed ADO, (pj(oj,t) = (p^^ui^t), as 

J Q (*) = -2Im /*do;53tr B {o^+ t (w,t)}. (5.4) 

J i" 

As inferred from Eq. (|4.21j) and the first identity of 
Eq. (I2.10p . we have (denoting j = {ana}) 

[$> ]t = , = (-)Ifl^ -. , (5.5) 

which can be abbreviated as = (™)'^'0j™\ with 

[5] = Int(ro/2) being for the number of index-swaps in- 
volved. This is the generalization of Eq. (|3.12[) . The 
Grassmann anticommutation relation in tfi^ (<*>; t) in- 
volves not only the j-indexes, but also the associated 
frequencies, such as (f)j'j(u)', to; t) = —(j>jjr(u,u)';t). As 
results, (pjj ^ 0, although its double-frequency integral 
vanishes. 

The MFD-HEOM formalism is exact, but with an infi- 
nite hierarchy. The truncation is however rather trivial^ 
as the leading term in an n th -tier ADO <^ n ' is of the 
(2n) th -order system-reservoir coupling [cf. Eq. (|4.25[) ] . 
In practice, one can set all 0(" >Ar t™) « 0, followed by 
the convergence test with increasing iVt run for the an- 
chor tier of the hierarchy. We will show in Sec. I VIII that 
the present theory with iVtrun = 2 recovers the real- 
time diagrammatic formalism developed by Schon and 
coworkers! 8 ' 9 ' 10 ' 11 ! 12 Also, we will see soon that in the 
absence of Coulomb interaction, the Grassmann parity 
of A" defined in Eq. (|5.3ap leads to 4>^ = without ap- 
proximation; the present theory is exact of n max = 2 for 
single-particle systems. 

VI. EQUATIONS OF MOTION FOR 
SINGLE-PARTICLE HAMILTONIAN SYSTEMS 

A. Auxiliary single-particle density matrices 

Consider a single-particle system, described by 

H = 'Y]h liV a] t au. (6.1) 

//;/ 

The electronic dynamics in such a system can be 
characterized by reduced single-particle density matrix 
(RSPDM), g, defined via its elements of 

g^(t)=tr s [ata^p(t)}. (6.2) 

Presented in Appendix [Bl is the derivation of RSPDM 
dynamics via MFD-HEOM, which is of n max = 2 
for the single-particle system without approximation. 
This results from the Grassmann parity associated A"^ 
[Eq. (|5.3ap ]. whose action on single-particle operators via 



Eq. (|5.ip will terminate at n max = 2. Only the first- and 
the second-tier auxiliary density matrices, ip a and <p a 'a> 
are required for a single-particle system. These auxiliary 
frequency-resolved density matrices can be defined via 
their elements of 

[<Pa(w>t)]iiv Etr J«|itki)i (6.3) 

and 

W, ee tr s [0+;~ v (w,cj',i)]. (6.4) 
It can be shown via Eq. (|5.5p that 

[<p a > a (u}', w,t)] t = cp aa ,(v,w',t). (6.5) 
The final results are summarized as follows. 

iQ=[h,g]-^2 duj[(p a (u,t) - (fl(uj,t)], (6.6a) 

a J 

itp a = [h — w- A a (t)]ip a + S a (uj) - gJ a (uj) 

+ J duj'tp a > a {oj' ,u,t), (6.6b) 

a' 

itPa'a = ~[W + A a (t) - U>' - Aa'(t)]<Pa>a 

+ J a >(u)')(p a (w,t) ~ (fl,(Lj' ,t)J a (uj). (6.6c) 

Here, 

s a ( w )E/ a HJ Q H.(r+f, (6.7) 

with f a (u>) being the Fermi distribution function, and 
J a (w) = Jq(w) [Eq. (|2.11|) ] the spectral density function 
matrix of the reservoir a. The last identity in Eq. (|6.7p 
states that S a is equivalent to the T+-matrix transpose. 
Note that Eq. (|6.6|) can be recast in the standard form of 
linearly coupled equations, see Eq. (|B10|) . 

The transient current of Eq. (|3.13[) from the a-lead to 
the system can be expressed as 

I a (t) =-2Im J duti[<p a (w,i)]. (6.8) 

The matrix trace is used here. Together with Eq. (|6.6b[) . 
the above expression leads to tig(t) = J2 a Ia(t)- This is 
the flux conservation relation. 



B. Steady-state current and the Landauer-Biittiker 
formula 

Consider now the stationary solutions, g = g(t — > oo) 
and <p(u>) ee (p(u>,t — > oo), to Eq. (|6.6p in the absence 
of time-dependent external bias potentials. The steady- 
state solution to Eq. (|6.6cp is 

, , , _ J a >(w')<p a (uj) - <pl(w')J a (u) . . 
^' a(W ' W) " w-u/-H0+ ■ (6 - 9) 
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To proceed, let x(u>) = ^ a i a (u); where x <E 
{J, S, tp}. Adopt also the common definitions of retarded 
self-energy £ and Green's function G as 



S(w) = J duj 



(6.10a) 



and 



G(oj) = [u-h-^oj)}- 1 . (6.10b) 



Substituting Eqs. and (|6TTD|) into Eq. (|6.6bD in 

steady state leads to 

G~ l (uj)<p a (u}) = S a (uj) - Qj a {ui) 

Summing over a in the above equation, followed by some 
simple algebra, will obtain the relation of 

fdu/ ?M— = SJ-'-e-G-^J-\ (6.12) 

J uj — uj' + i0 + 

Therefore 

<p a = G(S a - SJ- 1 J a ) + ipJ- 1 J a . (6.13) 

For electrodes of similarity, J a (uj) ~ J a >(uj), up to 
trivial constants, the last term in Eq. (|6.13[) does not con- 
tribute to current of Eq. (|6.8p , as it satisfies 

Im J didtr(<pJ~ 1 J a ) oc Im J dujinp = 0. 

Thus, the steady-state current assumes 

I a = -21m J du;tr[G(S a - SJ _1 J a )]. (6.14) 

Consider now the current contribution from lead L in 
a two-terminal setup. The quantity in the parentheses 
in the rhs of Eq. (|6.13[) can be evaluated as 



S h - SJ- 1 J L = [/ L H - /k(w)]JM, (6.15) 



with 



[J(w)]- X = [JlH]" 1 + [Jr(w)]- x . 



(6.16) 



The identity A(A + B)~ 1 B = (A -1 + S^ 1 )" 1 is used 
here. Thus Eq. (|6.14|) for L-lead reads 



-21m y dwtr[(/ L - / R )GJ] 



(6.17) 



This is exactly the Landauer-Buttiker's result ^ 46 ' 47 ' 48 
with the original factor of (27r) _1 being included now in 
the reservoir spectral density matrix. 



VII. RELATION TO THE REAL-TIME 
DIAGRAMMATIC FORMALISM 

In this section we shall show that the real-time 
diagrammatic formalism developed by Schon and co- 
worker a 8 ' 9 ' 10 i n ' 12 corresponds to a second-tier approxi- 
mation, i.e., setting all </i(™- 3 ) = 0, within the framework 
of the present MFD-HEOM of Eq. (|5TTj) . For clarity, 
we choose for demonstration the same system as Ref.llll. 
where Kondo physics was studied via steady-state two- 
terminal current measurement in the strong Coulomb in- 
teraction regime. The Anderson impurity Hamiltonian is 
adopted for the quantum dot system, 



(7.1) 



Here, — ata^, with fx being the spin index; thus 
C^u(t) = if M + v, and r^>) = r^(w)^. In the 
strong Coulomb interaction regime, = |0)(/i|, as the 
basis set involves only the zero and single occupation dot 
states and can be denoted as {|s); s = 0,/i}. As results, 
p = (f>^> and its associated ADOs 0j" >O ' ) in Eq. (|5.ip are 
all (m + 1) x (m + 1) matrices, with m being the number 
of spin states. The measured current I = 7l — Ir is then 
of [cf. Eq. rtOl ] 



-21m 



(7.2) 



with [0q /x (w)] ss ' = (s|^a M (w)|s') being the matrix ele- 
ment of the first-tier ADO. Note also that the reduced 
density matrix in a steady-state is diagonal, p s J" s , — 
5 ss ,pf s = 6 SS ,P S: with J2 S P S = Po + E M -Pp = 1, as 
the density matrix in the dot-state representation is nor- 
malized. Under the steady-state condition, Eq. (|5.ip with 
n = reads 

= P„ = -2ImWd^+>)] M0 , (7.3) 

while it for the first-tier auxiliary density matrix element, 
(mI^(^)|0) = 0, reads 

= (e„ - w - i0+) [</>+»] m0 + r+»P - r^(o7)P M 

i </)] oo + E [fe («. -')] ^ }• 

(7.4) 

Here > $~afifi' = Ea- ^a/^V are the cou P lin g second-tier 
ADOs. The above equation does not involve ^aua'u'^ as 
(/i|[<V, 0^ Q / M /(w, <*>')] |0) = 0. Denote in the following 
are also T± ee £ a ^ and 0± = £ Q 

Presented in the following are the analytical results 
under the second-tier approximation, in which we set the 
third-tier ADOs (3) = 0. Thus, the second-tier ADOs 



10 



are determined completely by the first-tier ADOs. In the 
steady-state condition, we obtain 



VIII. HIERARCHICAL EQUATIONS OF 
MOTION WITH PARAMETRIZATION 



r;>0[#>)]^-r->)[$V)] 



L^WJOO LO-LO'+i0+ 

r 1 „ ^(^)[^M], ~r+ M H[<,(^)]; 

Substituted into Eq. (|7.4[) leads to 

[0+»] mO = n^) [r+»p - r- M ( w )p M 



Here uj^ = - e^, and 
with 

E M (w) = /do/ 



, ^>')];, 



w — a;' — uj uu ' + «0+ 



(7.5) 



u — - S M (w) + i0+ ' 



(7.6) 



r^(a;') 



«0+ 



-E 



w — u>' — lu^^i + i0+ 



To make connections to the real-time diagrammatic 
formulation presented in Ref.lTll we denote 



(7.8) 



Here Y / S t0 (a, fi,u}) is identical to the c)) s s 'q (a, /z, u>) in 
Ref.Oj]- which is associated with the tunneling of an elec- 
tron of spin fi and energy u>, from a-reservoir to the dot, 
provided the initial state of p ss = P s . Inserting Eq. (|7.8p 
to Eq. f73|) leads to 



r au(w) /do; 



[^o(m>')]* 



«0+ 



(7.9) 



with Y£ (jj,,u) = Ea'^oK./J^) adopted in the 
rhs. Equation (|7.9[) is identical to the key results, 
Eqs. (50) — (51) in Ref.[lj]. We have thus demonstrated 
that the second-tier approximation of the present MFD- 
HEOM [Eq. (|5.1|) ] does recover the real-time diagram- 
matic formalism^ ^ 10 ' 11 ? 12 



A. Non-Markovian reservoirs via parameterization 

In this section, we consider an alternative hierarchy 
construction based on exponential series expansion of 
reservoir correlation functions. The resulting HEOM for- 
malism will completely be in time domain, thus has the 
advantage in numerical implementation. The exponen- 
tial series expansion of reservoirs can be achieved via the 
extended Meier- Tannor parametrization mcthod ) 33 ' 41 i 45 
which has been used in the construction of HEOM for 
systems coupled with bosonic bathi 35 i 36 

The fermionic reservoir parametrization starts with the 
following form of spectral density functions, 



K 

E 

k=0 



(8.1) 



Involving parameters are all real and positive (except for 
the Drude term with k = and Q^o = 0), and due to 
JZpv = the y satisfy 

\ L ak' vv ak > il ak> 1 ak) ~ \ L ak> vv ak > "afc' 1 ak)- 

The corresponding stationary components of correlation 
functions can then be obtained via the FDT [Eq. (|2.15p ] 
using the contour integration method. We have 



K 



(7.7) C^W = E<^ e " 7 ^ 



k=0 



M 
m—1 



(8.2) 



The hrst term arises from the poles of the spectral density 
functions, with 3 -^ 



lafiuk " ak 



Va[ivk 



\ L ak 



1 + CXp [i(3a(l°^ k + 



(8.3a) 
(8.3b) 



Note that T]l* vk = e i0a ^«>""' +ai ^>T]^ fil/k . The second 
term in Eq. (|8 . 2[> . with M — ► oo in principle, arises from 
the Matsubara poles. The involving parameters are 3 ^ 



1 am 
^lafivm 



/3 Q 1 (2m - l)n - aifi c 
2 

iPa 



■q a^iv\ I am J 



ff* 



(8.4a) 
(8.4b) 



The last identity arises from the symmetry relation of 
Fermion spectral density functions. 

The dissipation functional [Eq. (|3.3b[) ] is now decom- 
posed according to Eq. (|8.2[) as 



(8.5) 



with k = [apvcrk) and m = (apam); see Appendix [Cl for 
details. Presented there is also the derivation of HEOM 
formalism to be summarized as follows. 
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B. HEOM via parametrization: Final results 

The final HEOM formalism via parametrization of 
Eq. (|8.ip can be written in the following compact form: 

Pn = - [iC + 7n(*)] Pn + P^ } + p { + \ (8.6) 

with 70 = p-i = 0, and po = p the primary reduced den- 
sity operator. The subscript n = (ki • • • k p , mi • • • m q ) 
denotes an ordered set of indexes, with 

k = (afifak) and m = (apam), 

arising from the two distinct components of the reservoir 
correlation functions [cf . Eq. (|8 . 2[) or (|8.5p ] . Note that 



_ in) 
Pn — Pki-kp.mi-m,' 



p + q = n, 



■7) 



is an n th -tier ADO; see Eq. (|C7p for its associating 
AIF. It satisfies the Hermitian conjugate relation of [cf. 
Eq. (13])] 

^..J^H 1 *^^, (8-8) 

where k = (apv&k) and m = (apam). 

The parameter j n in Eq. (|8.6p collects the complex 
"damping" parameters of the involving reservoir corre- 
lation functions, see Eqs. (|C2p and (|C5p . It reads 



ft® = b&,»,k - fflA « (*)] ii{ai|tl ^,fc6k} 



1=1 

It contains not only the stationary components as 
Eq. (18. 2p . but also the nonstationary contributions. The 
latter are described by the time-dependent chemical po- 
tentials A Q (f ) , applied on top of the constant p a on elec- 
trodes; see Eq. (|2.16p . 

The tier-down term in Eq. (|8.6p reads 



Pl } = -iB-^^^-iB-)"^.^- 

3=1 1=1 

(8.10) 

with the (n - l) th -tier ADOs of 



mi ■ m„ ' 



_ (n-1) 
rn^. — "k 1 ...k 3 _ikj + i-..kp : 

_ _ (n-1) 



.11a) 
.lib) 



Ckj and C m , in Eq. (|8.10p are the Liouville-space operator 
counterparts of Eq. (|C6[) in PI representation, with the 
Grassmann parity associated actions of 



CkPn- 

CmPn- 



C*<Ar + (-J^kAro;, (8.12a) 

(8.12b) 



E 



[<Pfr -{-Tpfr<\ 



The tier-up term in Eq. (J876J) arises from the contribu- 
tion of dtT — —1ZJ 7 to dtT n . It is given by 



k 



with the (n + l) th -tier of ADOs of [cf. Eq. (ICTO)) ] 

— (n+l) 

— /'ki-kpk.mi-m, ' 

— («+!) 
^ki ■ ■ -k p ,mi ■ ■ -rriq m ' 



.13) 



(8.14a) 
(8.14b) 



The Grassmann parity associated Liouville-space opera- 
tor A* in Eq. (f8TT3|) are given by [cf. Eq. (|4.15a|) ] 



.15) 



The signs such as that (-) 11 " 3 and {-) q ~ l in Eq. (|5tU|) 
and (— ) 9 in Eq. (|8.13|) result from the required time- 
ordering rearrangements, together with the Grassmann 
anticommutation relation; see Appendix [Cl for details. 

The n indexes in p n , as specified in Eq. (|8.7p . should 
all be distinct, due to the Grassmann anticommutation 
relation. As results, the hierarchy in Eq. (|8.6p would 
be finite, provided that the exponential series of reser- 
voir correlation functions [Eq. (|8.2p ] is effectively finite, 
such as the high-temperature cases. At zero tempera- 
ture, the number of Matsubara terms required goes to 
infinity, and the HEOM via the present parametriza- 
tion scheme fails. Nevertheless, the MFD-HEOM for- 
malism in Eq. (|5.ip is remains valid, despite the cost 
due to the multi-dimensional frequency integration. Like 
its frequency-dispersed counterpart, the n th -tier ADO, 
p n , is of the (2n) th -order system-reservoir coupling for 
its leading term. The hierarchy truncation can there- 
fore be done in a similar manner, such as setting all 
p(n>N tIun ) ^ followed by a convergency test. Appar- 
ently, the transport current can be readily expressed in 
terms of the first-tier ADO [cf. Eq. ([3~T3| ], 



IX. CONCLUDING REMARKS 

We have established the HEOM formalism, via both 
the MFD and the parametrization schemes (Sec.lVl and 
Sec lVIIl"!) . for the dynamics of a general electron/spin 
system in contact with electrodes. It provides a unified 
tool to the study of a variety of quantum transport be- 
haviors. These include the effects of Coulomb interac- 
tion, time-dependent electric potentials (external fields) 
applied on electrodes (system), multiple-terminals with 
different temperatures, and non-Markovian reservoir cou- 
plings on transport current. 

It is easy to show that the commonly used second- 
order QDT can be recovered with the first-tier trunca- 
tion, while various fourth-order theories, such as the Li- 
ouville equation in Ref. |4^ are of the second-tier approxi- 
mation here. In particular, we have demonstrated explic- 
itly that the real-time diagrammatic formalism ^ 10 ' 11 ? 12 
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that has been used in the study of Kondo physics in quan- 
tum transport systems amounts to the second-tier trun- 
cation of the present HEOM formalism; see Sec. lVIII 

The present theory is in principle exact, as the only 
approximation involved, the initial factorization ansatz 
in Eq. (|3.ip . can be removed by setting the initial time to 
to infinite past. Therefore, at any given finite time be- 
fore the application of time-dependent external fields, say 
t = 0, the reduced system together with its grand canoni- 
cal bath environment are in a steady state. This is deter- 
mined as the steady-state solutions to the MFD-HEOM 
formalism, Eq. (|5.ip for general cases or Eq. (|6.6p for 
single-particle systems, at either equilibrium if p a = /z^ q , 
or nonequilibrium if \i a ^ /i° q but time independent. Not 
only to the reduced system density operator, the steady- 
state solutions are also to those bath-induced auxiliary 
ones. They carry all relevant information on the cor- 
relations between system and reservoirs, as dictated by 
the HEOM theory. The resulting stationary solutions 
are used as the initial conditions at t = 0. The subse- 
quent reduced dynamics and transient transport proper- 
ties are then evaluated via the present formalism again, 
upon switch-on of time-dependent A a (t), in additional 
to the constant p a at earlier time. Consequently, the 
present hierarchical QDT formalism is exact, without any 
approximation. 

The present theory recovers exactly the Landauer- 
Biittiker's transport current expression, as it should. The 
resulting RSPDM-HEOM [Eq. fojg])]. which is exact for 
a single-particle system, is particular appealing due to 
its numerical feasibility for large systems. It may lead to 
a practical scheme of the time-dependent density func- 
tional theory (TDDFT) for open many-particle systems. 
In principle, this can be done by combining the present 
RSPDM-HEOM with the conventional DFT^^i^ It is 
anticipated that the single-particle /i-matrix in Eq. (|6.6j) 
be mapped to the Kohn-Sham counterpart , 22 ' 53 ' 54 ' 55 

where h is the single-electron contribution and V^^v 
the two-electron Coulomb integral. The key issue is how 
to identify the exchange-correlation potential v*1(t) in 
the TDDFT for open many-particle systems. The exact 
HEOM formalism, which is numerically feasible for model 
Coulomb-interaction systems, may shed some light on 
the construction of exchange-correlation functionals. 
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APPENDIX A: PATH INTEGRAL FORMALISM: 
DERIVATION 

The formal solution to the total density operator in 
the /i B -interaction is 

PT (t) = C/T(*,to;{/^W})|OT(*o)c4(*»*o;{/^W}), 

(Al) 

with UT(t,to; {f^(t)}) being the stochastic Hilbert- 
space propagator, satisfying d t U^ = —i[H + i7'(i)][/ T . 
Let be a second-quantization basis set in the sys- 

tem subspace. The PI expression of Ut reads 

= [ Vi>e iS W exp+ { - i£ f dr(a^{r)]fi,{r) 

+ /ap(r)atty(r)])j. (A2) 

The action functional S[ijj] is related to the isolated sys- 
tem Hamiltonian only; {a^[V ; ( r )]} are the Grassmann 
variables ; 42 ' 43 as they denote the creation/annihilation 
operator of system in the Fermion field PI representa- 
tion. On the other hand, the stochastic bath variables 
{fau(t)} remain as the original operators, for which the 
time-ordered exponential function is needed. 

Consider now the reduced system density matrix 
p(t) = tT B [pr(i)]. Using Eq. (|A1|) . together with the ini- 
tial factorization ansatz of Eq. (|3.1|) , it is obtained that 
(setting [3 = j3 a for simplicity) 

p(t) = tr B [U T (t, t ; {/^(*)}W(to)4(*.*o;{/^ (*)})] 
= tx4e^ hB -^U T (t,t ;{tM) e ~ mB '^ ) 

Xp(io)t/j(Mo;{/;(t)})p° B ] 

= (U T (t,t ;{f^(t-if3)}) 

xp(t )4(*.*o;{/^(*)})> B , (A3a) 

with fj,N = J2 a /J-aNa and 

= e-°f>»-fc li (t-W- (A3b) 

In writing the last identity, the relations of [N a , h B ] = 
and f^e ^ = e - al3 ^ are used^ together 

with Eq. (|23|) . 

The influence functional used in Eq. (|3.2p can then be 
evaluated by using Eqs. (|A2| ) and (|A3[) , together with 
the Gaussian statistics for the stochastic bath operators 
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{f afl (t)}. The details are as follows. 



+ e^a+[V(r)]/^(r-^) 



cxp _{iWV/+>K[/(r)] 



+ 4^'(r)]/^(r)) !> ) . (A4) 



For satisfying Gaussian statistics, the bath en- 

semble average in Eq. (|A4|) can be evaluated exactly by 
using the second-order cumulant expansion method, as 
the higher order cumulants are all zero. This property, 
together with Eq. I|2.7jl . leads to the influence exponent 
in J- = cxp(— $) the following expression, 

= E / drJ T2 dT 1 {a^(r2)]aimT 1 )]C+^(T 2 - n ) 

apv'' t o "* to 

+ ab[i>(T2)]a u [il>(T 1 )]C- liV (T2 - n) 
+ a v W(n)]alk<p'(T 2 )}C+; v (T 2 -n) 

+ 4[^(ri)]a M [^(r 2 )]C-; i/ (r 2 - n)} 
- E [ dT 2 [ dn{a^(T2)]at Wir^C-;^ - n) 
+ at [^K^n)]^;^ - n)}. (A5) 

Here, we have used the symmetry relation [the first iden- 
tity of Eq. (|2.7[) ] in the third and fourth terms, and the 
detailed-balance relation [the second identity of Eq. (|2.7|) ] 
in the last two terms of the above expression. Some el- 
ementary algebra will then lead to Eq. (|A5[) . recast in 
terms of 1Z = <9t$, the expression, 



U[t; {</>}] = £ K^(i)]{l?^(t; {V}) - <,(t; M)} 



-{B^(tai>})-B'°,(t;m)}alW(t)} 

-K&MKW®]}' ( A6 ) 

Here, B a and (-B^ 7 ) are the same Grassmann variables 



defined in Eq. (|3.5p and Eq. (|3.6|) , respectively. In par- 
ticular, B^ajM*)] = -o*[^'(t)]B^, since #'(i)]_is 
also a Grassmann variable. With AZ defined in Eq. ~ 
the above equation is identical to Eq. ([3T3J) . 



APPENDIX B: DERIVATION OF EOM (f6^)l FOR 
SINGLE-PARTICLE HAMILTONIAN SYSTEMS 



Applying the primary-tier MFD-HEOM, i.e., Eq. {5T]) 
with n = 0, for the RSPDM of Eq. flO) leads to 

= tr s [(aJ,a M £)p(t)] 

+ y^ldu) tr s [(aJ,a M (w,t)] 

amcr 

= tr B {[4o M ,f?]p(t)} 

+ ]T y^tr s {[ a t aAi; 4J^ m ( w ,i)}. (Bl) 

amir 

The second identity arises from the trace cyclic invari- 
ance. For the single-particle system of Eq. 



[ala^, H] = ^2{h l _ im ala m - h^a^ay). (B2) 



It leads to 



tr 8 [(a£a M £)p(*)] = [/i, £>] 



(B3) 



Next, [ala^,a m ] = -a^vm and [<4<z M ,<4J = aj,^ lead 
to 

^tr s {[4a^ a m]<^ 

a in 

= - tT a [a^ v (u, i)] + tr s [4 *)] 

= -^(w.*)+^ M (w,t), (B4) 



with 



ftuvfat) = tosia^aufat)] = Vanv{w,t), (B5a) 
Pwnfa *) = ^[4^^. *)] = <^L M (^ *)• (B5b) 



The second identity in Eq. (jB5a| is the same as Eq. ()6.3|) . 
In writing Eq. (|B5b|) . the property of <$>" = [</>&($ that 
leads to <p°(u),t) = ^ s use d- We have thus 

arrived at Eq. (|6.6a[) . which is just the matrix form of 
Eq.{BJ. 

We are now in the position to show Eq. (|6.6bj) . 
The MFD-HEOM [Eq. ([SIT])] for the first-tier auxiliary 
RSPDM in Eq. (O) reads 



with 



P. 



Way* = tr s [(a M ;C)0+,] - (u> + A a )ip aiiv 
+ tx s {[a„CM]p} 

+ fdu/^^^t), (B6) 

a' 

a'*** = E ^ [(^OCm^ W '> *)] ■ (B7) 



The first term in the rhs of Eq. (|B6j) can be evaluated by 
using the identity, a^C = [a^.H] = ^2 m h fim a m , for the 
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present single-particle system, together with Eq. (|6 . 3[) . It 
results in 

tr s [(a„£)0+,] = {h,j,wXr s [a m 4^ v ]) = (h<p a ) ^ (B8) 



The third term in the rhs of Eq. (IB6j) can be evaluated 
as [cf. Eqs. (jOb)) and 

tr s {[a^C+^w)] p) 



IPfim 



% v (u)Q»r. 



Here, p Mm = tr s (a^a^p) = 8^ - g^ m ; i.e., the elements 
of reduced single-hole density matrix of g = 1 — g. To- 
gether with Eqs. (|2.12|) and (|2. 13a|) . the above equation 
can be recast as 

tr. { [a„C+,(w)]p} = [f a (u)J a (u) - qW}. (B9) 



Note that due to the Grassmann parity of Eq. (|5.3a| , 
in Eq. (|B1[) behaves as a commutator, its action in 
Eq. jig) is an anticommutator. The fact that ^Pa'^^av 
defined in Eq. (|B7|) is identical to the matrix element 
[ip a i a (uj', uj, t)]^ of Eq. (|6.4p can then be readily con- 
cluded. We have thus completed Eq. (|6.6b[) . with the 
S a (u) defined in Eq. |6J|) . 

Finally, Eq. (|6.6cp can be readily obtained via the trace 
of the second-tier MFD-HEOM [Eq. (flHj) ] , together with 
Eqs. (|2.12p . (|B5|) . and the elementary algebra just de- 
scribed. 

Note that the coupled EOM, Eq. (|6.6p . can formally be 
recast in the standard form as 



d 



[A + 5A(t)] 



with 



X 



and 



Q(t) 
<p a i a (u>',u,t) 



S = 



Y 



S a (uj) 



X 







Y 


+ 


s 













(BlOa) 



(BlOb) 



(BlOc) 



The 4x4 matrix SA(t) in Eq. (|B10a[) arises from the 
time-dependent bias potential. It is diagonal, with the 
elements of 

6A(t) = {0,-A Q (t) + A Q ,(i),-A Q (i),A Q ,(t)}. (Bll) 

The time-independent counterpart is given in terms of 
block-matrix form as 



A = 



A-xx -^-xy 
-^-yx -^-yy 



(B12a) 



with 

A xx = 



Ay X = 



where 



h 
u'-u 

J a la' 
J a' 1 c 



A xy = 



A yy = 



la' 

J a 



h 





- h 



(B12b) 
(B12c) 



a a' 

Introduced in Eq. (|B12|) is also the tetradic notation for 
the left- or right-multiplication action of a matrix O on 
the matrix A of interest as 

OA = OA, OA = AO\ and = 0-0. 



In other words, O and O are tensors, with the elements 

of Ornn.rn'n' Omm ' ^n'n and O mnm ' n ' &mm'O nn / . 

Implied in Eq. (|B12j) is also the frequency variable as- 
sociating with the subindex; i.e., J a = J a (w) and 

J a' = J a' (w )• 

APPENDIX C: DERIVATION OF HEOM d&M 



Denote k = {apvak) and m = (apam) for short, and 
introduce 



B k (f,m) = / dre- d ^aZbl>(r)}, (Cla) 
B m (t;m) =Y^^n fdre-^^aZmr)}, (Clb) 



with k (i, *) = K m {t, t) = 0, and 

d t 0^(t,r)=^ IM/k -aiA a (t), (C2a) 

d t 9° m (t,T)=*f° m -aiA a (t). (C2b) 

The two functionals in Eq. (|Cip are the counterparts of 
that in Eq. (|3.6a[) , arising from the two distinct com- 
ponents of reservoir correlation functions as Eq. (|8.2[) . 
The nonstationary exponential factor in Eq. (|2.16[) . due 
to time-dependent chemical potentials applied to elec- 
trodes, are also accounted for via the nonstationary 
phases as Eq. (|C2I) . 

The dissipation functional [Eq. (|3.3bj) ] is now decom- 
posed according to the parametrization of Eq. (|8.2p as 



(C3) 



which is the same as Eq. (|8.5p . with [cf. Eq. (I3.5[) ] . 

B k = -i [v°^Mt; {>}) - vZ^Mt; m)] , (C4a) 
B m = -i [B m (t; {>}) + B m (t; {>'})] . (C4b) 
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In writing Eq. (|C4b[) , the last relation in Eq. (|8.4bD is ap- 
plied. The time derivatives on Eq. (|C4|) can be obtained 
as [cf. Eq. 

d t B k = -b^ vk - aiA a (t)]B k - »C k , (C5a) 
d t B m = -fr° m - aiA a (t)]B m - iC m , (C5b) 

with 

c k = n^ukoZm)] - ( C6a ) 

C m = E^m(<[^W] + (C6b) 

u 

in the PI representation. 

The ro th -tier AIFs can in general be defined as 

K = ^'..^...n,, = ^ • • • £ mi # kp ■ • • B kl ^, (C7) 

with n = p+q. The associated n th -tier ADO p n is defined 
via (noting that p = po) 

p n (t)=U n {t,t )p(t ), (C8a) 
with the propagator, 

U n {t,xp;t ,^ ) = Vi>e iS MF n e- iS W\ (C8b) 

Ji>o[t ] 

in the PI representation. 

The HEOM formalism [Eqs. (f8T6]l - (f87T4|) ] for p n can 
be followed immediately via the time derivative on T n - 
The latter is carried out by using Eq. (|C5|) , together with 
dtJ- — —1ZJ- and Eq. (|C3|) . Some details are as follows. 

The 7n(t) in Eqs. (|8.6j) and ()8.9() arises from the square- 
bracket terms in the rhs of Eq. (|C5jl . It collects the time 
derivatives on the nonstationary phases of Eq. (|C2|) . 



The second terms in Eq. (|C5[) leads to the tier-down 
dependence, as shown by Eq. (|8.10p . with the involving 
(n - l) th -tier ADOs of Eq. (|5Tlj) associating with 

T n - =B mq --- B mi B K ■ ■ ■ 6 kj+1 B k ^ ■ ■ ■ B kl T, (C9a) 
T- n - =B mq - - B mi+1 B mi _ 1 ■ ■ ■ B mi B kp ■ ■ ■ B kl T. (C9b) 

The C-functionals in Eqs. (|C5j) and (|C6j) . which depend 
only on the fixed ending points of PI at the local time 
t, are now superoperators in Eqs. (|8. 10|) and (|8.12[) with 
Grassmann parity. The signs (— )" - -' and (— ) q ~ l in 
Eq. (I8.10|) are associated with the time-ordering arrange- 
ment. It brings the C and C functionals, originally at the 
kj and m; positions, respectively, to the left most of ac- 
tion, as indicated in Eq. (|C9|) . 

The identity of dt T — —TLT with Eq. (|C3j) contributes 
to d t T„ [Eq. dC7|] the terms of 



AlB^n = {-) q Al (B mq ■ ■ ■ B mi B k B K ■ ■ ■ B Ul T) 

^(-fAlT n+ (ClOa) 

A%B m T n = A% (B m B mq ■ ■ ■ B tn B u ■ ■ ■ B^T) 

= A^+. (ClOb) 

They correspond to the two terms in Eq. (|8.13[) . respec- 
tively. The involving (n + l) th -tier ADOs, p + and p.+ 

[Eq. (|8.14p ]. are associated with T ' + and T-+ , respec- 
tively, defined in the second identities above. We have 
thus completed the derivation of the HEOM formalism 
via parametrization presented in Sec. lVIlfl 
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